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Abstract. A new method for analyzing point patterns 
produced by the evolution of gravitational clustering is 
presented. The method is taken from the study of molec- 
ular liquids, where it has been introduced for making a sta- 
tistical description of anisotropic distributions. The statis- 
tical approach is based on the spherical harmonic expan- 
sion of angular correlations. A general introduction to the 
method is given; a theoretical analysis shows that it is par- 
tially connected with previous harmonic analyses applied 
to galaxy catalogs. The effectiveness of the statistical anal- 
ysis in quantifying clustering morphology is illustrated by 
applying the statistical estimators to point distributions 
produced by an ensemble of cosmological N— body sim- 
ulations with a CDM spectrum. The results demonstrate 
that the statistical method is able to detect anisotropics in 
non-random patterns, with different scales being probed 
according to the expansion coefficients. 
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1. Introduction 

According to the standard picture the observed struc- 
tures in the Universe have arisen via gravitational in- 
stability from the time evolution of the initial matter 
density fluctuation field S m (x). At very early epochs the 
field is assum ed to be a random Gaussian process (e.g. 
Peebles 1980). In the proposed scenarios the formation 



of structures proceeds hierarchically, with galaxies form- 
ing first and larger structures following later from clus- 
tering of galaxies. The present galaxy distribution shows 
that clustering is ranging from small groups up to clus- 
ters and superclusters. A striking feature of the observed 
large-scale galaxy distribution, as revealed by extended (~ 
100 - 200Mpc) reds h ift surveys flGeiler fc Huchra 1989 ; 
da Costa et al. 1994; Landy et al. 1996; Broadhurst et 
ah 1990; (Saunders et al. 1991a ), is that galaxies are 
connected in a web of sheet-like structures, with 
voids between them a nd filaments at the intersections 
flGeiler fc Huchra 1989| ). The size of these structures can 
be as large as 50 ~ lOQMpc. These features revealed in 



the galaxy distribution have been formed through gravi- 
tational clustering and are connected to the initial power 
spectrum P(k) =< \S m (k)\ 2 >; thus the observed galaxy 
clustering can in principle be used to put constraints on 
the allowed cosmological models, for which the fundamen- 
tal parameters determine the shape of P(k). The galaxy 
distribution today is highly non-linear and N— body meth- 
ods have been used to sample the initial phase space dis- 
tribution and to study the time evolution of the cluster- 
ing growth ( Efstathiou et al. 1985| ; Bertschinger & Gelb 
1991). These methods lead to a final particle distribu- 
tion which should be a representative sample of the ex- 
pected galaxy distribution for the cosmological model un- 
der study. The use of N— body codes implicitly makes the 
assumption that galaxies trace the matter distribution. In 
order to compare the observed galaxy clustering with the 
results of numerical simulations one has to introduce a 
statistical descriptor. The method used to investigate the 
statistical distribution of galaxies is therefore very impor- 
tant because it should be used to discriminate between 
different models, furthermore the method must also pro- 
vide a mathematical description for the rich variety of 
clustering morphology as seen in the galaxy distribution. 

The first method to be introduced as a tool for study- 
ing galaxy clustering was the 2-point correlation function 
£(r) ([Totsuji fc Kihara 1969fc peebles fc Hauser 19740. 



Analyses of galaxy catalogs show that £(r) has a power- 
law shape £ oc r -7 , with 7 = 1.77 and is equal to 
unity for r Q = (5 ± l)hT x Mpc Q flPccblcs fc Hauser 1974 ; 
Fisher at al. 1994). The 2-point correlation function has 



been widely used to constrain cosmological models from 
results of numerical simulations (Jenkins et al. 199S and 
references cited therein). However the 2-point function 
does not provide useful information about the rich vari- 
ety of structures which characterize the morphology of the 
galaxy distributions. The reason is that the power spec- 
trum P(k) ( the Fourier transform of £ ) does not describe 
the correlation of the phases that the 5k 's develop during 
the clustering evolution. 



1 Here Ho = WOhKmsec 1 Mpc 1 is the present value of the 
Hubble constant 
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For a complete description of the clustering one has 
in principle to consider correlation functions £jv of higher 
order. The measurement of these correlations is imprac- 
tical however when N > 5 because of the noise in- 
troduced by the finite number of galaxies in the sam- 
ple (Szapudi & Colombi 1996). For this reason differ- 
ent statistical tools have been introduced for study- 
ing the clustering of the galaxy distribution. These 
methods are : the void probabili ty distribution Pp(V ) 
(White 1979; |Vogeley et al. 1994 iGhigna e t al. 1997|) ; 
moments of counts in cells ( [Saunders et al. 19911 ); prob- 
ability distribution of co unt in cells ( Bouchet et al. 1993 ; 
Ucda fe Yokoyama 19961) . 

While these methods still make use of the higher 
order moments of the clustering distribution, alterna- 
tive approaches have been introduced which are more 
geometrical and implicitly contain information from 
the N— point correlations at all orders. The first 
method to be introduced was the percolation statis- 
tic QZerdovich 1982|; [Klypin fc Shandarin 1993|), other 



methods include the topological genus-density threshold 



( pott, Dickinson fe Melott 1986| ; Ryden et al. 1989), min- 
imal spanning trees (Barrow, Bhavsar & Sonoda 1985); 
Minkowski functionals (Mecke, Buchert & Wagner 1994; 
Kcrscher et al. 1997| ), graph-theory ( |Ucda fc Itoh 1999| ) 
or a global descriptor (J-function) based on the nearest- 
neighbor distribution ( Kerscher at al. 1999|) . In order to 
quantify the geometrical features of the galaxy clus- 
tering several statistics have been explicitly developed 
to single out in a quantitative way the filaments and 
sheets of the cosmic network. These statistics are based 
on the moments of the mass distribution and can be 
used to detect filamen ts ( Fry 1986 ; Dave et al. 1999 ; 
Babul fc Starkman 1992 ) , o r the shape of the cluster - 
ing ( |Luo fc Vishniac 199$ |Robinson fc Albrecht 1996] ). 
The shape statistic can also be constructed us- 



ing Minkowski functionals ( ^chmalzing, fc Buchert 1997 
Bharadwaj et al. 2000"|). 



The existence of this variety of approaches is related 
to the fact that there is not any single well defined theo- 
retical method which can be used to analyze the shape of 
the clustering network. Furthermore a statistical method 
must also be able to produce a robust statistical measure, 
in order to be applied to galaxy catalogs. These require- 
ments are naturally satisfied by the 2-point function £, 
but, as outlined above, this function it is not indicated for 
detecting filaments. 

In this paper, I propose an alternative method for an- 
alyzing the clustering morphology which is a generaliza- 
tion of the 2-point function £ and is based on the spher- 
ical harmonic analysis. The statistical method has been 
applied in molecular fluid dynamic simulations for de- 



metallic glasses ( 


Steinhardt, Nelson & Ronchetti 1983; 


Wang & Stroud 1991 


). The method is applied here for il- 



lustrative purposes to point distributions obtained from 
cosmological N— body simulations. 

The paper is organized as follows: In Sect. 2, I present 
the method. Sect. 3 describes the ensemble of CDM cos- 
mological simulations used to test the effectiveness of the 
method. In Sect. 4 the statistic is applied to the point sets 
produced by the simulations, the main results are dis- 
cussed and the conclusions are summarized. 



2. THE METHOD 

2.1. Spherical harmonic analysis 

For describing galaxy clustering evolution in the non- 
linear regime, use is made of different statistical methods, 
based on the galaxy positions, which have been outlined 
in the Introduction. The approach introduced here makes 
use of the positions of the point distribution within a cer- 
tain distance from a randomly chosen one. 

The method is drawn from molecular dynamics simu- 
lations, where it has been introduced for studying orienta- 
tional order of supercooled liquids and metallic glasses. A 
general review can be found in Haile & Gray (1980) and 
Mc Donald (1986), here I will follow the notation of Wang 
& Stroud (1991). 

Let us consider a system of N p particles. The i—th par- 
ticle has coordinates r^, in an arbitrary reference frame. 
For a specified cutoff radius R c all the particles such that 
|r, — Vj\ < R c are neighbors of i. 

The line joining i to one of the j is termed a bond. 

The angular coordinates of the vector A.ji = Tj — Ti 



are 9j , <fij and the quantity 



Qlm(ri) = £ Yi m (0j,(/)j), 



(1) 



is the coefficient of the spherical harmonic expansion of the 
angular density of the bonds associated with the particle 
i. 

In Eq. [I], and hereafter, summation is understood over 
all particles j of the distribution such that \ri — Tj\ < R c - 

The coefficients Qi m (ri) are defined as the bond- 
orientational order parameters (Wang & Stroud 1991) 



and they can be drastically changed by a rotation of the 
reference systems. A natural quantity to consider , which 
is rotation invariant, is 



\ 7U— — 1 



(2) 



Using the addition theorem for the spherical harmon- 
ics, the expression simplifies to 



(3) 
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where Pi is the Legendre polynomial, 

ljk = cos(6 jk ) = Aji • Afci/dA^HAfcil) 

is the angle between two bonds and the summations in 
Eq. |^ are then independent of the chosen frame. 

The expression for Qi(ri) can be further simplified if 
one considers that jjj — 1, furthermore it is more conve- 
nient to redefine Qi{ri) dividing by the mean number of 
neighbors < N% >= N of the particle i : 



(4) 



j k>j 



For a random distribution the summation terms in 
the square root are negligible for Ni 3> 1, thus Qi{ri) oc 
for Poisson noise. 

For the whole system an order parameter < Qi > can 
be defined by averaging over all of the N p particles: 



x n p 

< Qi >=jf^Qi{n)=Qi- 

P i=l 



(5) 



Thus the clustering distribution can be studied by eval- 
uating the Qi at a specified cutoff radius. Going to higher 
coefficients (I 3> 1) means probing smaller scales in the 
bond distribution. 

However the Qi do not exhibit spatial information, and 
a more useful quantity is the auto-correlation Gi(r) of the 
coefficients Qi(ri). The function Gi is defined as follows: 
for all of the M p pairs (i, k), such that \ri —rk\ = r ± Ar, 
where Ar is the thickness of the radial bin, then G;(r) is 
the sum over all of these pairs 



1 



M„ 



47T x - 

21 + 1 ^ 



Qim(ri)Qi m (ri +r). (6) 



This equation can be greatly simplified: let j be the 
set of neighbors of the particle i and p that of the particle 



k, which satisfy \rj — 
the summation becomes 



< R c and |r„ — rd < R c . Then 



QLi r i)Qlm(ri+r) 



E 



m=l 

E E^( f 

m=-l j 

The final expression for Gi is: 

g <m^EEEEw p )> 



(7) 



with Tj p being the angle between Aji and A p/ t . The sum- 
mation over the pairs is X^iEfc! w hh the sum over the 
particles k only for those particles with jr^ — = r±Ar. 



Note that now it is not possible to simplify the equa- 
tion as was done for Eq. ^, because one is considering the 
correlation between different order parameters. 

The information about the degree of angular correla- 
tion between bonds contained in the Gi (r) is more explicit 
if the Gi themselves are divided by the zero-order corre- 
lation function 

Go(0 = — E E Oo(n)Qo(r< + r). (8) 

k 



Mr, 



The functions Gi(r) can be made dependent also on the 
direction of r by considering all the pairs with orientation 
of the relative separation r in the solid angle range £l r ± 
ASl r . These functions Gi(r) are of scarce practical use 
because of the noise level which is introduced when applied 
to a finite sample. 

2.2. Theory 

The equations that have been obtained for Qi and Gi(r) 
are defined through a point set distribution; from a theo- 
retical point of view it is possible to relate these functions 
to the power spectrum of the density fluctuations. Spher- 
ical harmonic analysis has been widely applied in cosmol- 



ogy, both to the galaxy distribution (Bcharf et al. 1992 
[Scharf ct al. 1993 ; Ballingcr, Heavens & Taylor 1995) and 



to the study of cosmic background radiation maps 



( |Hu, Sugiyama fc Silk 19971 ; [Bartlctt 1999[ ). I will follow 
here the notation of Scharf et al. (1992), who have an- 
alyzed the angular distribution of the IRAS redshift sur- 
vey using spherical harmonics. For the generic distribution 
considered in Sect. 2.1 let us define n(r) to be the number 
density and n to be the average one. It is also useful to 
introduce <i>(r) as a generic selection function, here $ = n. 

Then the coefficients Q" l (r a ), r a being the center po- 
sition, are given by 

QT( r o) = f $(O^17"0V)dV r = r + r', 
J n 

where the integral is taken over a sphere of radius R c . The 
number density is 

n(r) — n[l + S(r)] — n[\ + S(r + r )], 

here S(r) is the density fluctuation. Thus the equation for 
QT( r °) becomes 



<5>[l + 5]Y, m d 3 r = 

$Y, m d 3 r' + ( <S>{r')5{r + r' )Y, m d 3 r' , 



Q?{ro) 



where the first term on the rhs is zero for I > 1. The den- 
sity fluctuation S(r) can be expanded into Fourier modes, 
and the final expression for Q" l (r D ) is 
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The e %kr is now decomposed into spherical functions, so 
that the ensemble average gives 



<fr(x)x 2 ji(kx)dx 



where x and y are the particle separations relative to the 
origin and £ is the reduced three-point correlation func- 
tion. Integration over solid angles removes for / > 1 the 
first three terms on the rhs of the equation. Of the remain- 
ing terms, the first ( £(|sc - y\ = J2k |<5fc| 2 e lfe ( x ~ !/ )), gives 



(10) 



here the ji(x)s are the spherical Bessel functions. The 
summation over k can be translated into an integral and 
JlY^fdn = 1, so that 



2 f°° f^ c 
< |Q™| 2 > = V- / P(k)k 2 dk / §{x)x 2 ji{kx)d. 
^ Jo [Jo 

+Af, 



(11) 



where P(k) =< \Sk\ 2 >, V is a normalizing volume and 
A/" is the shot noise term TV = J $(a;)2; 2 cfcc. 
The coefficients C; are defined as 



C, 2 = 



1 



A/"(fi/47r)(2/ + l) 



m |2 
Z I > 



(12) 



and assuming statistical isotropy they can be related to 
Eq.(||) by Qi ~ Ci/N^ 2 . For an uncorrelated distribution 
C; = 1, furthermore the integral J [| Rc 3>(x)x 2 ji(kx)dx = 
^i(kR c ) defines a bidimensional window function in the 
parameter space (R c ,l). Increasing the value of I means 
probing smaller angular scales, the value of the cutoff ra- 
dius R c sets the size of the spectrum contribution to th e 
integral (see, for example, Fig. 4 of Scharf et al. 199 . 
The statistical method introduced in Sect. HA, and origi- 
nally proposed in a different field, is then mathematically 
equivalent to the spherical harmonic analysis which has 



been applied to galaxy surveys (Scharf et al. 1993). The 
order parameter Qi is just related to the angular power 
spectrum estimator C; ( Scharf et al. 1992] ). 

However Eq. (|l0|) is valid as long as the position of 
the observer located at r can be considered as random 
with respect to the bond distribution inside the sphere of 
radius R c . This is a good approximation for values of the 
cut-off radius R c much larger than the 2-point clustering 
length, a condition which is valid for the analyzed angular 
catalogs. In Sect. 4 the statistical method is tested using 
the point distributions obtained from a set of cosmologi- 
cal simulations, and the coefficients Q\ are calculated for 
different parameters and cut-off radii R c . Therefore one 
has to consider in principle also the three-point correla- 
tion terms between the particle located at Tq and the 
other particles in the bond sphere. The ensemble average 
for \Q\ n \ 2 is then (Peebles 1980, sect. 36B): 



< \QT\ 2 > = / *dx J $dy[i + e(a) + H(y) + £(l* - v\) 
+a^yAx-y\)]Yi m {n x )Yr n {n y ), (is) 



Eq. fllpp, while the other takes into account for the trian- 
gle configuration the correlations between the three parti- 
cles. The integral can be evaluated in terms of the Fourier 



transform of Q (Fry 1984), the bispectrum B(k±, k 2l k 3 ) : 



J <Pdx J <Pdy((x,y,\x~y\)Y l m (n x )Yr i (n y ) = 

k 1 k 2 k 3 

el (k lXl +k 2X2+ k 3X3 ) = B {k 1 = -(fc 2 + fe 3 ), fe 2 , fc 3 ) 



k 2 k 3 



$dx J $dyy ; m (fi x )r i * m (f7 y )e l ( fc2X+fc33/ ) 



(47r) 2 H'^ / / B(fci,fc 2 ,*!3)*l(*i^)*l(*2^c) 

k 2 k 3 

Yr(n kl )Yr^k 2 ) , (w) 



where the constraint ^ ki = is required by homo- 
geneity, x = x 2 — xi and y = x 3 — x\. In second 
order perturbation theory, for Gaussian initial condi- 
tions, the bispectrum can be expressed as ( Fry 1984 ; 



Verde, Heavens & Matarrese 2000) 



B(k!,k2,k 3 ) = K(k 1 ,k 2 )P{kl)P(k2) + cyc. 



(15) 



here eye. means cyclic terms. The term AT(fci,fc2) 

is weakly dependent on the cosmology and can 

be written as (Fry 1984; G oroff et al. 1986; 
[Verde, Heavens fc Matarrese 2000 ): 

K(ki,k 2 ) = A Q + A lC os(9 12 ) + A 2 cos 2 {9 12 ), (16) 

where the coefficients Ai depend on the assumed biasing 
expansion for the local density field and 9 12 is the angle 
between fci and k 2 . The integral (Q) can then be evalu- 
ated for a particular cosmological model. The dependence 
of B on the triangle shape in the fc-space shows that the 
integral is non zero only for I < 2. The integral (|lj) is a 
particular case of the integrals considered by Verde, Heav- 
ens and Matarrese (2000). These authors have expanded 
the projected galaxy density in spherical harmonics and 
studied the 3-point function of the expansion coefficients, 
which is a quantity directly related to the bispectrum. 
The evaluation of the above integral is rather complicated 
(see, e.g., Eq. (28) of Verde, Heavens and Matarrese 2000) 
and will not be considered here. In Sect. 4 it will be seen 
that the corrections ( |l4| ) to Eq. (10) must be rather small 
already for R c ~ 2tq. 

This theoretical analysis can also be used to derive 
an expression for the quantity G/(r), which contains spa- 
tial information about the bond angular distribution. The 
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evaluation of Gi(r) now involves correlation terms up to 
the fourth order. In analogy with Eq. ([l3]) one has to com- 
pute for Gi(r) the ensemble average 

Gi(r) =<Qy>(r )Qt m (r + r)>= 
J <$>dx J <f>dy[£(\r +y-x\) + ((x, \r + y-x\,\r + y\) 
+C(\r +y-x\,y,\r- x\) + £(|r + y\)£(\x - r\) + £(r) 
£(|r + y-x)\) + r)(x, r,y,\r + y\,\r - x\,\r + y - x\)\ 
Yr(n a )Y? m (Q v ), (17) 

where r\ is the reduced four-point correlation function, and 
only the terms which are non-zero after the integration 
over the solid angles are shown in the integrals; the sep- 
arations of the four-particle configuration of Eq. (Q) are 
defined as: r = r k - n, y = r p - r k , x = rj - r*. 
The integrals in (Jjjj) can be evaluated using the same 
method adopted for the integral (|l4]). As previously dis- 
cussed, from the results of Sect. 4 for < \Q] n \ 2 > , the third 
and fourth-order terms are expected to be negligible. I will 
consider here only the leading order term £(|r + y — x\) 
and the ensemble average for Gi(r) is then 

<Q?{r )Qr{r + r) >= 



,— ikr 



|y,"WI 2 [j^ Jl (kx)x 2 dx\ 2 (A^) 2 \8 k \ 2 . (18) 



The summation over k becomes an integral and if 



-ikr 



is now decomposed as 



e -ikr = 4 n Y,t P 3 P (kr)Y p * q (n r )Y pq (n k ), 



(19) 



then the triple integral for the spherical harmonics Y™ in 
Eq.pl) is 



jY m Y^Y lm d£l k = / " Y m {-) m Yi_ m Yi m dQ 
{-) 



(21+iy (2 P +i) 

4tt 



1 / 2 ( I I p 





I I p 

in —m q 



( h h h \ 
The symbols are the Wigner 3—7 coef- 

V mi m-i rtij, J 

ficients. Because of their properties, the summation over 

the q terms in Eq. ( |l8| ) is non-zero only for q = m — m = 0, 

furthermore the p summation admits only terms with 

21 + p — even, < p < 21. Finally the expression for 

Gi(r) becomes : 

47T 



Gi(r) 



21 



1 J2<QT(r )Qr(ro + r)> = 



V I dkk 1 
o 



<frx 2 ji(kx)dx 



P(k) 



(27t) 3 

p(even)— '' 



m=—l 



i i p\( i i p 

000/lm-mO 



(20) 



In the above equation the function Gi still depends on 
the orientation of r through the m terms. However, ow- 
ing to statistical isotropy, this dependence should not 



E m (-) r 



be present. In fact, it can be shown that (Gangui 1995) 
I I p 
m —m 
20|) simplifies to 



(21 + l)( 1 /2)(_)/^ Qi and Eq . 



Gi(r) 



(2ttY 



V 



dkP(k)k 2 



f Q " $(x)x 2 ji(kx)dx j (kr 



(21) 



From this function one can recover the standard 
definition for £(r) taking for I = the limit R c — * 
0, J^ c <&(x)x 2 jo(kx)dx — > 1/47T. The statistical descrip- 
tor of clustering defined by Eqs. (||) and (0), when applied 
to a spatial point process, should be able to describe in 
a quantitative way the patterns generated by the cluster- 
ing distribution. As outlined in the Introduction different 
methods have been employed to this end; the most impor- 
tant advantage of the order parameters introduced here is 
that they can can be used to study spatial patterns in 
a well-defined way by varying the angular scale / or the 
depth R c of the bond spheres. In order to test the clus- 
tering features described by the order parameter Qi and 
correlation function Gi , I have then analyzed point dis- 
tributions generated by a set of cosmological TV-body sim- 
ulations with CDM spectra. The simulations are evolved 
in time starting from an early epoch when linear theory 
applies. The clustering evolution gives rise to a variety of 
spatial structures and the statistical descriptors are then 
applied, in order to analyze the produced patterns, to nu- 
merical outputs selected at different epochs. The results 
obtained are discussed in Section 4. 



3. THE SIMULATIONS 

The point distributions used for the statistical analysis 
have been constructed from a set of purely gravitational 
cosmological TV-body simulations. For each of these simu- 
lations the clustering evolution of N p — 10 6 particles is fol- 
lowed in time in a box of comoving size L = 200h~ 1 Mpc. 
The gravitational potential is solved using a P 3 M code 
(|Efstathiou ct al. 1985|), with N m = 256 3 grid points. The 



interparticle force has a linear shape cloud profile with a 
comoving softening parameter e — 0.220Mpc. 

The initial conditions for the particle positions and 
velocities have been generated at the initial time ti 
according to the standard Zcl'dovich approximation 



(Efstathiou et al. 1985). The power spectrum P(k) is that 
of a standard CDM model ( Bardeen et al. 198(f ) with 
Cl m = 1 and h — 0.5. The Fourier components 6 k are 
Gaussian distributed, with random phases and variance 
P(k). An ensemble of five different integrations has been 
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carried out, each with a different random realization for 8k- 
The simulations are integrated in time with the expansion 
factor a(t) being used as the integration variable and with 
a(ti) — 1. Records of particle positions have been taken at 
a(t) = 3.0,3.8,4.5,6,7.9. The power spectrum is normal- 
ized so that the rms mass fluctuation in a R — Rh^Mpc 
sphere takes the value <t$ — 0.4 at a(t) = 3.0. 

In order to measure the order parameters associated 
with peaks above a given height of the underlying den- 
sity field 8(x) I have used the peak-background formalism 



0.3 



(White et al. 1987). For each particle i the peak number 
Hi is calculated at a(ti) from the initial density field; ni 
is the number of peaks with height 8 S > va s , here 8 8 
is the field 8(x) Gaussian smoothed with R s = l.lMpc 
and al is the variance of 8 S . The field is constrained to 
take the value Vb&b when smoothed on a scale i?t, > R s 



flWhite et al. 1987| ; |Park 199l|) The particle i is then iden- 
tified as a 'peak' if rii > p, with p being a uniform random 
number between and 1. 

The total number of particles associated with the 
peaks of a given height is ~ n i- Three different thresh- 
old levels have been chosen: v = 0.5, 1.3, 2.5 which corre- 
spond, respectively, to ~ 250,000, 165,000, 33,000 
in the simulation box. Hereafter the notation {x} s means 
the subset of the 10 6 simulation particles selected to iden- 
tify the peaks with threshold v = s. At different epochs 
a(t) the order parameters have been computed by apply- 
ing the statistical estimators to the point distribution pro- 
duced by {x} s . The clustering of the mass distribution 
produced by all of the simulation particles has been also 
considered, in this case 8 > s = — 1. 



4. RESULTS AND CONCLUSIONS 

The order parameters Qi are shown in Figs. 0, | & | 
for different values of the expansion factor a(t), thresh- 
old levels and cut-off radius R c . In Fig. [j] the time evo- 
lution of Qi is plotted as a function of I for different 
epochs. Several generic features can be seen in the plot: 
the coefficients Qi grow with time as the clustering devel- 
ops and decay with increasing / as smaller angular scales 
are probed. This follows directly form the shape of P(k). 
For a random test distribution with Ni ~ 10 3 neighbors 
Qi ~ const ~ l/y/N, and the estimator (||) is clearly 
able to detect non-random distributions with a high level 
of significance. The scatter in the numerical ensemble can 
be judged from the size of the error bars, for the sake of 
clarity they have been shown in each plot only for a single 
case. For the plot of Fig. |] the point distribution is that 
obtained from the {^}-i set. In this case the computa- 
tion of the Qi can be greatly reduced by using a random 
sub-sample of the particle set. This has been done also for 
other sets {x} s , according to required computational load 
needed to evaluate the coefficients Qi. 

The results obtained have been found to be quite ro- 
bust to changes in the size of the random sub-sample, usu- 
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Fig. 1. The coefficients Qi are plotted as functions of 
I for different expansion factors a(t). The estimator 
has been applied to the whole point distribution of the 
simulation particles. Here and in the following plots the 
cut-off radius R c is in units of h~ 1 Mpc. Error bars show 
the scatter within the five simulation ensemble. The con- 
tinuous line with a nearly constant value for the Qi refers 
to a random test distribution. 



ally an accurate evaluation of Qi can be obtained already 
with a ~ 20% random sub-sample of the particle set. How 
the Qi change when particle populations with different 
clustering properties are selected is shown in Fig. ^. For 
a fixed a(t) and R c the coefficients Qi have been evalu- 
ated for particle subsets associated with different thresh- 
olds. For any I the coefficients Qi get larger when higher 
thresholds v are selected. The choice of higher v 1 s corre- 
sponds to a more clustered distribution, and this is clearly 
detected by the coefficients Qi. It is worth stressing that 
the trend is obtained according to the adopted definition 
(^) for Qi(ri). The angular power spectrum coefficients C\ 
scale approximately as cx Ni, thus the enhanced clustering 
of particle subsets associated with higher v' s is masked for 
the Ci by the drop in N, with respect lower thresholds. 
The black squares in Fig. |^ are the rescaled linear the- 
ory coefficients (|l2|). In this case the comparison has been 
made for the population 8 > — 1, in order to avoid errors 
introduced by a linear bias approximation for Eq. ([l2|), if 
subsets {x} s corresponding to higher thresholds v would 
have been selected. 
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Fig. 2. Here the estimators Qi have been computed at 
a fixed time and cut-off radius R c from the point distri- 
butions generated by {x} s , with s = —1,0.5, 1.3,2.5. The 
black squares are the rescaled linear theory coefficients C/ , 
computed from the CDM power spectrum at a(t) = 6 and 
for R c = 20. The power spectrum has been normalized 
as in the simulations, so that the points compare directly 
with the 5 = — 1 Qi coefficients. The black triangles refer 
to the same theoretical coefficients, but computed in this 
case using the non-linear power spectrum of the simula- 
tion. 

There is a good agreement with the numerical experi- 
ment up to I ~ 5; beyond this angular scale ( ~ tt/5) non- 
linearity effects make the Qi depart from the linearized 
coefficients. For the adopted normalization the size of the 
non-linear scale is about Rnl — 5 — 8h~ 1 Alpc, at the 
epoch shown in the plot, which is in rough agreement with 
R C A9. When the power spectrum of the simulation parti- 
cles is used to compute the theoretical coefficients there is 
an excellent agreement with the corresponding measured 
Qi down to the smallest angular scale (black triangles of 
Fig. |^) . For / < 2 the small differences show that the cor- 
rections ( pl| ) to Eq. jlC|), for the considered cosmological 
model are neglig ible already for R c > 2Qh~ 1 Mpc. 

The coefficients Qi are plotted for different cut-off ra- 
dius i? c , in Fig. |3| for a fixed time and a given set {x} s . 
There is a general tendency for the Qi to get smaller at a 
fixed I when R c is increased. This is expected because the 
bond distribution approaches isotropy as R c gets higher. 
This trend is confirmed for all values of I and choices of 



Fig. 3. The coefficients Qi are shown at a(t) = 6 for 
four different cut-off radii R c = 10, 20, 40, 60 ( in h~ x Mpc 
units). The estimators have been applied, as in Fig. |, to 
the whole particle distribution. 



different distributions {x} s . Note also how the scattering 
in ensemble is reduced (R c — 40h~ 1 Mpc) with respect to 
the smaller values of R c shown in the previous plots. The 
analysis performed so far shows that the order parame- 
ters Qi can be used as reliable statistical indicators of the 
global clustering properties of a point distribution. 

From the histograms of Fig. || an interesting result is 
that the quadrupole coefficient Qi of the bond angular 
distribution is larger than the dipole coefficient Q\. This 
is valid with a high significance only for values of the cut- 
off radii R c ,$ 2ro, close the value of the 2-point cluster- 
ing length. This excess of local anisotropies, compared to 
what is found for R c £ 20h~ 1 Mpc, suggests that a sig- 
nificant amount of substructure is present in the particle 
distribution of dark matter halos. Jing (2000) has ana- 
lyzed the halo dark matter density profiles from a set of 
cosmological N— body simulations. For those halos which 
are substructure rich and less virialized, he found large 
deviations from the analytical profile of Navarro, Frenk 
& White (1995). The above results therefore suggest that 
the proposed statistical method can be profitably used to 
correlate in a quantitative way the local halo anisotropies 
with other halo properties, when investigating their evo- 
lution in different cosmological models. 
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Fig. 4. The bond correlation functions Gi(r), divided by 
Go, are shown at a(i) = 6 and R c = 40 for the distribution 
obtained from {a;} s =o.5- The different histograms are for 
different values of I, only Gi up to I = 4 are plotted. 
The line with a value for Gi close to zero is for a random 
distribution. 



To analyze the clustering shape, more interesting re- 
sults are obtained if the estimator is applied to point 
distributions generated by different {x} s . The total oper- 
ation count of applying Eq. ([?]) to a given point distribu- 
tion scales as oc N*, thus making it prohibitive to evaluate 
Gi(r) for the whole set {x} s of the simulation box. There- 
fore one has to resort to compute Gi taking as centers i 
in Eq. (R) a random subsample of {x} s . The summations 
over the particles k at a distance r from i and the neigh- 
bor particles of the chosen centers ■ J2 P ) are instead 
run over the whole set {x} Sl although a random dilution 
of the neighbor sets was performed in many cases in order 
to reduce the total computational cost. For a given parti- 
cle i taken as a center the mean number iV, of neighbors j 
within a distance R c from i ranges from ~ 10 2 to ~ 10 4 , 
according to the chosen cut-off radius R c , threshold level 
v and expansion factor a(t). 

However the use of a random subset requires some care 
with respect to the previous procedures used to compute 
the coefficients Qi. It has been found that the function 
Gi(r) is subject, in some cases, to large sample-to-sample 
variations when different random subsets are chosen from 
a given {x} s to evaluate the corresponding G/. This shows 
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Fig. 5. As in Fig. fibut in this case Gi (r) have been eval- 
uated summing in (Q) over all the neighbors of the center 
particles; in the previous figure G;(r) were calculated tak- 
ing random sub-sample of the neighbor sets. 



that Gi(r) is a sensitive measure of the clustering patterns 
generated by a point distribution. 

One of the consequences of this result is that the calcu- 
lated Gi (r) analyze in this case the clustering morphology 
of the point distribution as it is measured according to 
the selected random sub-sample, rather than that of the 
whole {x} s . In order to consistently compare G/(r) for 
different choices of several parameters, the same random 
subset has been chosen at different times and for different 
R c . Furthermore G;(r) for point distributions obtained 
for different thresholds have been calculated with random 
subsamples subject to the constraint that for a given {x} s 
the random subset was part of the random sample cor- 
responding to a lower threshold, and so on down to the 
lowest level with s = — 1 . The application of these con- 
straints leaves about 8, 000 particles in the simulation vol- 
ume, which have been taken as center particles when the 
G;(r) were calculated. This procedure allows a consistent 
analysis to be made of the clustering morphologies de- 
tected by the Gi (r) at different epochs and different values 
of the cut-off parameter R c , angular scale I and threshold 
v. The results obtained are shown in the following plots. 
Hereafter G\{r) is the estimator (R), applied to point dis- 
tributions, divided by Gq{t). Because of the large number 
of cases spanned in the parameter space, only a few of 
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Fig. 6. The functions Gi(r) plotted here have been ob- 
tained from different point distributions and the same val- 
ues of /, R c and a(t) (given in the panel). The distributions 
are those given by the simulation particle subsets {x} s , 
with s = 0.5, 1.3, 2.5, and {a;} s= _i. 



them have been plotted, the main features seen are quite 
general. 

In Fig. [| Gi (r) is given as a function of r for different 
values of I, up to I = 4. The values for the other pa- 
rameters are given in the panel. The correlation between 
angular distributions for different bonds decreases with in- 
creasing r and as smaller scales are probed, but with some 
exceptions. In this case Gi, for 1 = 1, has a significant 
correlation above the noise extending for a wide range of 
spatial distances. In other cases a similar trend has been 
found also for I — 2. Gi(r) is a measure of the degree of co- 
herence in the clustering patterns. For the example shown 
in the figure the result is a web of large-scale structures, in 
the examined point distribution, which is detected by the 
estimator. These features are detected at the lowest order 
only (I < 2), since their distribution in the bond spheres is 
limited to large angles. For r — > the correlation strength 
increases because the two neighbor spheres overlap and 
partially trace the same structures. The continuum limit 
( ^l| ) is reached in the large sample regime Ni — > oo. For 
a finite number of points Gi (r) measures the anisotropies 
of the distribution which samples the density field. 

A theoretical estimate of how Gi(r) are affected by 
noise terms is difficult, therefore an indirect check has been 



Fig. 7. The Gi(r) are shown for different values of the 
cut-off radius: R c = 10, 20, 30, AOh^Mpc. The values of 
the other parameters are given in the panel. 



performed for the examples showed in the plots evaluating 
Gi (r) in correspondence of a random distribution. For each 
of the considered cases this was obtained by considering 
tha same set of M p pairs used in Eq. (m) to evaluate Gi (r) , 
but with a random distribution of N r — 10 3 neighbors 
within each of the two spheres of radius R c centered at 
the pair coordinates. The results show that Gi ~ for a 
random distribution of neighbors, so that the measured 
Gi(r) are well above the noise. 

For the example showed in Fig. || Gi (r) have been com- 
puted using random sub-samples of the N ~ 10 4 neigh- 
bors of the center particles i. In this case the fraction / 
of selected neighbors was / » 1/4. Finite sample effects 
on the measured G;(r) can be estimated by repeating the 
evaluation of Gi(r) for the case of Fig. ^|, but this time 
running the summations in (Q) over all the neighbors par- 
ticles. The results are showed in Fig. ^ a comparison with 
the previous figure shows that the Gi (r) profiles are quali- 
tatively similar, suggesting that measurements of the clus- 
tering pattern by Gi(r) are not affected by a dilution of 
the neighbor sets. Empirically it has been found that for 
the examples showed here this is valid for fN XL 10 2 . 

Gi(r) is also sensitive to different degrees of clustering, 
as shown in Fig. |^. In this case the G; have been computed 
for different {x} s with the other parameters being kept 
constant (see the panel). 
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Fig. 8. The correlation functions Gi(r) are plotted here 
for I = 1,R C — 40, a(t) — 6 and point distributions ob- 
tained from {ir} s= 2.5. Different histograms refer to {x} s 
from simulations with different random realizations of the 
ensemble. 



It is important to stress that this result was obtained 
because of the way in which the random subsets were cho- 
sen; if different subsets had been used for different {x} s , 
then sample-to-sample variations would have dominated 
over the differences in Gi for different thresholds. A generic 
feature of Gi is a decrease in the signal as R c is increased 
(Fig. [?]). This is expected, since in going to higher cut-off 
radii, the point distribution in the bond sphere approaches 
isotropy. 

As previously stressed Gi (r) has been found to be very 
sensitive to sample variations, the size of the fluctuations 
can be checked in Fig. ||, where G/(r) for a specified set 
of parameters has been computed for the different {x} s 
of the numerical ensemble. Because of the computational 
cost of evaluating Gi (r) I have calculated these functions 
for the whole simulation ensemble only for a(t) = 6. The 
plot shows clearly how fluctuations in the ensemble, in 
this case, are of the same order as the variations in Gi 
shown in the previous plots when the input parameters 
were changed. 

In one simulation large-scale correlations are clearly 
detected for r ^ 60/i -1 Afpc. These results are in accor- 
dance with those of Doroshkevich, Fong & Makarova 
(1998), who have analyzed the evolution of large-scale 
structure in standard CDM TV-body simulations. Their 



Fig. 9. The same as in Fig. 
tained from {x} s 



but for distributions ob- 



=0.5- 



analysis of sheets and filaments in the particle distribu- 
tion is based upon the 'core sampling' method (Buryak, 
Doroshkevich & Fong 1994). A proper comparison is dif- 
ficult because of the different statistical methods em- 
ployed. The results obtained by Doroshkevich, Fong & 
Makarova (1998) reveal a richness of structures on scales 
of ~ 10 — \QQh~ 1 Mpc, in agreement with the correlation 
range measured here by the functions Gi. 

In Fig. ^ the functions Gi (r) are shown, for different 
realizations, for the same example as in Fig. but in this 
case for {x} s= o.5 instead of {x} s= 2.5- The scatter in the en- 
semble is clearly reduced; this is a consequence of the lower 
threshold considered. For a point distribution approach- 
ing that of the whole mass in the simulation volume then 
Gi(r) has a larger contribution from particles which are 
not part of filaments or of other large-scale structures. The 
signal is thus reduced as also are the variations between 
different random realizations. This is valid for cut-off radii 
much larger than the non-linear scale. If R c = \Qh~ 1 Mpc 
then the functions Gi would have shown no sign of correla- 
tions, with strong fluctuations around zero independently 
of the separation distance r. The scattering in the ensem- 
ble is also reduced when smaller angular scales are probed 
(Fig.®. 

To summarize, the observed results have shown that 
the statistical method defined by the functions Qi and 
Gi(r) can be used to analyze the clustering morphology 
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Fig. 10. As in Fig. §, but for I = 6. 



cause of the tensor window function of the angular mask 
arising from partial sky coverage. The required formal- 
ism has already been applied in the literature to angular 
catalogs (Scharf et al. 1992; Scharf et al. 1993) and its ap- 
plication to the Qi coefficient is straightforward. Further 
complications arising from the radial selection function of 
the catalog and redshift space distortions have also been 
considered (Scharf et al. 1993; Ballinger, Heavens & Taylor 
1995). For the functions Gi(r) the analysis is much more 
cumbersome and is beyond the scope of this paper, which 
is intended to introduce the method and to illustrate its 
features by applications to numerical simulations. 

An issue which has not been considered is whether this 
statistical tool can successfully be used to check the consis- 
tency of clustering data with cosmic models of structure 
formation. The analysis performed reveals (Fig. |J) that 
any method which quantifies in a sensible way the pres- 
ence of filaments or planar structures in the cosmic web 
is also prone to cosmic variance. When different scales or 
levels of bias are considered, then the fluctuations are re- 
duced (Figs. H & [h]). The results obtained suggest that 
an application of the estimator G; (r) to real data will re- 
quire different scales to be probed, in order to give useful 
constraints on cosmic theories. 
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produced by gravitational clustering in a quantitative way. 
The method was originally devised in a completely differ- 
ent field for studying anisotropic structures and its ap- 
plication to cosmological clustering shows that, to some 
extent, it overlaps with previous analyses. The usefulness 
of the method has been tested by applying it to a set 
of cosmological N-body simulations with a CDM power 
spectrum. An application of the estimators to selected nu- 
merical outputs shows that the statistic is clearly able to 
discriminate between particle populations with different 
degrees of clustering. Several shape statistics have been in- 
troduced to quantify the presence of filaments or sheet-like 
structures in the clustering network. The function Gi{r) 
used here describes anisotropies in the clustering distribu- 
tion by measuring the degree of correlation between the 
angular densities as seen from two different observers sepa- 
rated by r. Gi(r) can then be considered a statistical mea- 
sure of clustering patterns, with different scales probed by 
varying the input parameters I and R c . This provides an 
alternative to other methods proposed so far to quantify 
the clustering morphology. 

With large redshift surveys becoming available in the 
next few years, the proposed statistical method appears 
as a promising tool for analyzing patterns in the galaxy 
distribution. In practice real galaxy catalogs have an in- 
complete sky coverage which implies, for the observed har- 
monic coefficients, non-zero mean values and a coupling to 
different modes of the all-sky coefficients. This occurs be- 
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